Simulation scripts/Coital Frequency/Network.Estimation.R

In EvoNetHIV/Goodreau_et_al_Behavior_-_SPVL: Evonet: Integrated Bio-Social Models for HIV Epidemiology

#script to 1) estimate and 2) save to file a "network" object to use in evonet
#this step can not be done in hyak and must be done on linux and then object transferred
# to hyak
#
# NOTE! Very important!!!!
# "network" related parameters used to estimate network here need to be same on hyak run
# These "network" parameters include (if different from default)
# initital population size, target stats, relationship duration, role ,
# generic attribute, nw_form_terms, nw_coef_form, nw_constraints

#--------------------------------------------------------------
library(evonet)

#--------------------------------------------------------------
# change warnings to errors and open browser on error
options(warn=2) # turn warnings into error
options(error=browser) # go into debug mode on error

# change error setting back to default if desired
#options(error=NULL)
#options(warn=1)

#--------------------------------------------------------------
#Load default parameters

primary_parameters  <- input_parameters_primary()
cd4_data            <- input_parameters_cd4_data()

#--- combine individual parameters into single list
evoparams <- c(primary_parameters, cd4_data)
#--------------------------------------------------------------
target_stats_vec <- c(0,250,500,750,1000)
coital_freq_vec  <- c(0.1,0.15,0.2,0.25)

modname_vec=c("conc0","conc5","conc10","conc15","conc20")
modname_vec2=c("msa0.1","msa0.15","msa0.2","msa0.25")

for(ii in 1:length(target_stats_vec)){
  for(jj in 1:length(coital_freq_vec)){
evoparams$initial_pop       = 5000
evoparams$initial_infected  = 500
evoparams$birth_model       = "poisson_birth_numbers"
evoparams$poisson_birth_lambda = 0.685
evoparams$trans_RR_age         = 1.0
evoparams$target_stats         = 5000*0.7/2
evoparams$relation_dur  			      	= 200
evoparams$trans_RR_insertive_anal_msm = 2.9  
evoparams$trans_RR_receptive_anal_msm = 17.3
evoparams$transmission_model  	      = "hughes"
evoparams$nw_form_terms = "~edges + concurrent + offset(nodematch('role', diff=TRUE, keep=1:2))"  

evoparams$target_stats   <- c(1750, target_stats_vec[ii])
evoparams$mean_sex_acts_day <- coital_freq_vec[jj]

modname=modname_vec[ii]
modname2=modname_vec2[jj]
#add parameters that are functions of other input parameters
evoparams  <- input_parameters_derived(evoparams)

#convert raw parameter list into EpiModel object
evoparams <- do.call(EpiModel::param.net,evoparams)


#--------------------------------------------------------------
# initialize network

nw <- setup_initialize_network(evoparams)

#--------------------------------------------------------------

#run qaqc on input parameters
input_parameters_qaqc(evoparams)

#--------------------------------------------------------------

#estimate initial network (create argument list, then call fxn)
netest_arg_list <- list(
  nw            =  nw,
  formation     =  as.formula(evoparams$nw_form_terms),
  target.stats  =  evoparams$target_stats,
  coef.form     =  evoparams$nw_coef_form,
  constraints   =  as.formula(evoparams$nw_constraints),
  verbose       =  FALSE,
  coef.diss     =  dissolution_coefs( dissolution =  ~offset(edges),
                                      duration    =  evoparams$relation_dur,
                                      d.rate      =  3e-05) )

estimated_nw <- do.call(EpiModel::netest, netest_arg_list)

#--------------------------------------------------------------

save(estimated_nw,
     file = paste("evo_nw_",modname,modname2,".RDATA",sep=""))
remove(estimated_nw)
  }
}
# change error setting back to default
options(error=NULL)
options(warn=1)